Assignment for EE5731 AY2023/2024 Semester 2

Exploring Edge and Corner Detection from Ultra-dark Images and Finding Epipolar Lines from Stereo Image Pairs (100)

Background

While understanding of vision information through computer vision techniques have been developed (and still developing) over the past decades, current techniques may not be applied directly to real-world applications due to various constraints. One of the concerns is that current techniques are mainly developed assuming the input data are shot under a non-challenging environment, with adequate illumination and contrast. This leads to the observable fragility of current techniques, which may not be capable to generalize well to adverse environments, including dark environments with low illumination. Take security surveillance as an example: obtaining crucial visual information such as edges and corners could play a vital role in detecting objects for anomaly detection. However, anomalies are more common at nighttime and in dark environments, yet current computer vision approaches are largely obscured by darkness, and are unable to extract key visual information effectively. It is therefore highly desirable to explore how we could obtain visual information given the current computer vision approaches in ultra-dark images.

Requirements

In this assignment, you are to complete several tasks that explores ways to obtain edges and corners from ultra-dark images, and to find epipolar lines from stereo dark image pairs.

Before doing the project, please read the instructions carefully (failure to follow them may result in unnecessary penalty):

1. Implement your codes with the code block in jupyter notebook, you may add in however much code blocks you require. TODO code blocks are added as placeholders. DO NOT modify the instructional markdown blocks;
2. Make sure your codes clean, easily readable (add meaningful comments if needed, comments are markable);
3. Make sure that the submitted notebooks have been run and the cell outputs are clearly visible.
4. Write your answers in the given markdown cells (with [Fill in xxxx]), keep your answers clear and concise;
5. Your submission should be the exported PDF file of your jupyter notebook. It is easier for you to "Export to HTML" first, then Print to PDF. Please print your PDF strictly in Portrait mode. Name your PDF file strictly to "MATRICULATION NUMBER_YOUR NAME_CA1.pdf"
6. Incorrect submission format will result in a direct 20 points (out of 100) deduction.
7. Do submit your project on Canvas before the deadline: 5:59 pm (SGT), 4 October, 2024;
8. Policy on late submission: the deadline is a strict one, so please prepare and plan early and carefully. Any late submission will be deducted 10 points (out of 100) for every 24 hours.
9. This is an individual project, do NOT share your solutions with others, we have zero tolerance for plagiarism.

Before Task 1

• For the following task, you are to use Pic1 as the ultra-dark image (denote as Dark-1). The corresponding long-exposure normal brightness image is GT1 (denote as Bright-1).
• All raw images are given in RGB format (colored images), but you may also convert to grayscale images for your convenience. Show the step and result of grayscale conversion first if you are to convert to grayscale for your tasks.
• IMPORTANT! You may use any function of OpenCV or its equivalence for basic operations (e.g., loading images, matrix computation, etc.), but strictly NOT the direct functions for each individual task/step (e.g., cv2.Canny or its equivalence for Canny edge detection, and cv2.equalizeHist for histogram equalization). Using such functions would consider the code to be erroneous.

Task 1: Canny edge detection and Harris corener detection in ultra-dark images (60%)

In this task, you will need to implement the Canny edge detection algorithm and try to improve its performance with image enhancement methods. You are to discuss how the characteristics of ultra-dark images affect the performance of Canny edge detection, and how different image enhancement methods would improve/affect the performance of Canny edge detection on ultra-dark images. You are to follow the following steps (requirements):

Detailed Steps/Requirements for Task 1:

1. Compute the image histograms of both Dark-1 and Bright-1, then discuss your observations, what is the characteristics of an ultra-dark image and how it is observed through the image histograms. (10%)
2. Perform Canny edge detection on both Dark-1 and Bright-1. For at least one of the images, show the output of every single key step of Canny edge detection (e.g., after image smoothing, after computing partial derivatives, after applying threshold, etc.) as well as the final output. Observe and compare the performances of Canny edge detection on Dark-1 and Bright-1 and discuss why you would observe your result. (10%)
3. Perform Harris corner detector detection on both Dark-1 and Bright-1. For at least one of the images, show the output of every single key step of Harris corner detector detection as well as the final output. Observe and compare the performances of Harris corner detector detection on Dark-1 and Bright-1 and discuss why you would observe your result. (10%)
4. Implement histogram equalization (HE) from scratch to Dark-1. Output the results of HE (denote as Dark-1-HE) and discuss your observation: how HE change/improve the sampled images. Hypothesize how will the output of Canny edge detection and Harris corner detection be like for Dark-1-HE. (10%)
5. Apply your implementation of Canny edge detection and Harris corner detection onto your HE-enhanced sampled images and demonstrate the final output. Do the final outputs fit your hypothesis? If yes, rationalize your hypothsis. If not, describe and explain the differences. Also discuss how are the results compared to that of Bright-1? (10%)
6. Lastly, choose a image enhancement method you prefer. State what image enhancement method has been chosen. Implement it with appropriate comments and output the results (denote as Dark-1-Self). Apply Canny edge detection and Harris corner detection onto Dark-1-Self and display the results. Discuss the differences between the Dark-1-Self against Dark-1-HE. Further, observe, compare, and rationalize the difference between the edges and corners detected between Dark-1, Dark-1-HE, Dark-1-Self, and Bright-1. (10%)

• Note for Step 6: you may use open-source codes or direct functions in OpenCV or equivalent for the chosen image enhancement method.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 1 IS HERE                           ##
#############################################################################

from PIL import Image
import matplotlib.pyplot as plt

def calculate_brightness_histogram(image_path):
    # load the image
    image = Image.open(image_path)
    grayscale_image = image.convert("L")
    width, height = grayscale_image.size
    # calc the histogram
    histogram = [0] * 256
    for y in range(height):
        for x in range(width):
            brightness = grayscale_image.getpixel((x, y))
            histogram[brightness] += 1
    return histogram

image_path_dark = 'Pic1.jpg'
image_path_bright = 'GT1.jpg'
image_pathes = [('dark', image_path_dark), ('bright',image_path_bright)]

histogram1 = calculate_brightness_histogram(image_pathes[0][1])
histogram2 = calculate_brightness_histogram(image_pathes[1][1])

# draw histogram
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
fig.suptitle('brightness histogram')

axs[0].bar(range(256), histogram1, color='blue')
axs[0].set_title('Dark Image')
axs[0].set_xlabel('brightness value')
axs[0].set_ylabel('pixel count')


axs[1].bar(range(256), histogram2, color='red')
axs[1].set_title('Bright Image')
axs[1].set_xlabel('brightness value')

plt.show()

[Fill in your discussion to Step 1 here]

Dark Image Observation

It seems like the pixels of dark image mostly concentrate in the low brightness value area, and the crest of dark image pixels also lays on the extremely low brightness value area.

Bright Image Observation

While the the pixels of Bright image distrubute more evenly and has crest laying on mid-low brightness area, which has lower height.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 2 IS HERE                           ##
#############################################################################
import cv2
import numpy as np
import matplotlib.pyplot as plt

# First create a canny detection function
def canny_edge_detection(image_path, image_name, low_threshold=35, high_threshold=100, show_steps=False, show=True):
    """
    Arguments:
    - image_path: path to the image file.
    - low_threshold: low threshold for the hysteresis procedure.
    - high_threshold: high threshold for the hysteresis procedure.
    - show_steps: whether to show intermediate steps of the algorithm.

    return:
    - edges: the final edge image.
    """
    # read the image
    image = cv2.imread(image_path)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    if show_steps:
        plt.imshow(gray, cmap='gray')
        plt.title(image_name + ' ' + 'original gray image')
        plt.axis('off')
        plt.show()
    
    # gaussian blur
    blurred = cv2.GaussianBlur(gray, (5, 5), 1.4)
    if show_steps:
        plt.imshow(blurred, cmap='gray')
        plt.title(image_name + ' ' + 'image after Gaussian blur')
        plt.axis('off')
        plt.show()
    
    # gradient calculation
    grad_x = cv2.Sobel(blurred, cv2.CV_64F, 1, 0, ksize=3)
    grad_y = cv2.Sobel(blurred, cv2.CV_64F, 0, 1, ksize=3)
    # gradient magnitude and angle
    magnitude = np.sqrt(grad_x**2 + grad_y**2)
    angle = np.arctan2(grad_y, grad_x) * (180 / np.pi)
    angle[angle < 0] += 180  # modify angle to [0, 180]

    if show_steps:
        plt.imshow(magnitude, cmap='gray')
        plt.title(image_name + ' ' + 'gradiant magnitude')
        plt.axis('off')
        plt.show()
    
    # Non-Maximum Suppression
    def non_max_suppression(magnitude, angle):
        H, W = magnitude.shape
        output = np.zeros((H, W), dtype=np.float32)
        angle = angle % 180

        for i in range(1, H-1):
            for j in range(1, W-1):
                q = 255
                r = 255

                # 0 degree
                if (0 <= angle[i,j] < 22.5) or (157.5 <= angle[i,j] <= 180):
                    q = magnitude[i, j+1]
                    r = magnitude[i, j-1]
                # 45 degree
                elif (22.5 <= angle[i,j] < 67.5):
                    q = magnitude[i+1, j-1]
                    r = magnitude[i-1, j+1]
                # 90 degree
                elif (67.5 <= angle[i,j] < 112.5):
                    q = magnitude[i+1, j]
                    r = magnitude[i-1, j]
                # 135 degree
                elif (112.5 <= angle[i,j] < 157.5):
                    q = magnitude[i-1, j-1]
                    r = magnitude[i+1, j+1]

                if (magnitude[i,j] >= q) and (magnitude[i,j] >= r):
                    output[i,j] = magnitude[i,j]
                else:
                    output[i,j] = 0
        return output

    suppressed = non_max_suppression(magnitude, angle)
    if show_steps:
        plt.imshow(suppressed, cmap='gray')
        plt.title(image_name + ' ' + 'Image after non-maximum suppression')
        plt.axis('off')
        plt.show()
    
    # Double Threshold
    strong = 255
    weak = 50

    output = np.zeros_like(suppressed, dtype=np.uint8)
    strong_i, strong_j = np.where(suppressed >= high_threshold)
    zeros_i, zeros_j = np.where(suppressed < low_threshold)
    weak_i, weak_j = np.where((suppressed >= low_threshold) & (suppressed < high_threshold))

    output[strong_i, strong_j] = strong
    output[weak_i, weak_j] = weak

    if show_steps:
        plt.imshow(output, cmap='gray')
        plt.title(image_name + ' ' + 'Image after double thresholding')
        plt.axis('off')
        plt.show()
    
    # Hysteresis Thresholding
    def hysteresis(img):
        H, W = img.shape
        for i in range(1, H-1):
            for j in range(1, W-1):
                if img[i,j] == weak:
                    if ((img[i+1, j-1] == strong) or (img[i+1, j] == strong) or (img[i+1, j+1] == strong)
                        or (img[i, j-1] == strong) or (img[i, j+1] == strong)
                        or (img[i-1, j-1] == strong) or (img[i-1, j] == strong) or (img[i-1, j+1] == strong)):
                        img[i,j] = strong
                    else:
                        img[i,j] = 0
        return img

    edges = hysteresis(output)
    
    # Final result
    if show_steps:
        plt.imshow(edges, cmap='gray')
        plt.title(image_name + ' ' + 'Canny Edge Detection Image')
        plt.axis('off')
        plt.show()
    
    if not show_steps and show:
        plt.imshow(edges, cmap='gray')
        plt.title(image_name + ' ' + 'Canny Edge Detection Image')
        plt.axis('off')
        plt.show()
    
    return edges

edge_dark_1 = canny_edge_detection(image_path_dark, "Dark-1", low_threshold=10, high_threshold=40, show_steps=True)
edge_bright_1 = canny_edge_detection(image_path_bright, "Bright-1", low_threshold=35, high_threshold=120, show_steps=True)

[Fill in your discussion to Step 2 here]

Observation

Generally, It seems that the canny edge detection performs better on bright picture than on dark picture. The edge picture of bright-1 has more clear edge on bright area while the dark-1 edge picture has nearly no edge in dark area and less clear edge on bright area.

Note: The double threshhold is seperately adjusted for both dark-1 and bright-1, however it's still hard to improve the dark-1's performance under only canny detection.

Possible Reasons

• The pixels of gray picture of dark-1 mostly lay on the area of extremely dark area and are easy to be filtered by low threshold in the stage of 'Double Thresholding'.
• The weak edge pixels are too much due to the extreme low light, and make it hard to find nearby strong edge in the stage of 'Hysteresis Thresholding'.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 3 IS HERE                           ##
#############################################################################
import cv2
import numpy as np
import matplotlib.pyplot as plt

def harris_corner_detection(image_path: str, name: str, k=0.04, threshold=0.01, show_steps=False, show = True):
    """
    realizer Harris corner detection and show each intermediate step according to the show_steps parameter.

    arguments:
    - image_path: Path to the image file.
    - k: Harris corner sensitive coefficient [0.04, 0.06].
    - threshold: Threshold for corner detection.
    - show_steps: Whether to show intermediate steps of the algorithm.

    return:
    - corners: Image with detected corners marked in red.
    """
    # Read the image and convert it to grayscale
    image = cv2.imread(image_path)
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    gray = np.float32(gray)
    if show_steps:
        plt.imshow(gray, cmap='gray')
        plt.title(name + ' origin gray image')
        plt.axis('off')
        plt.show()
    
    # compute gradient Ix and Iy
    Ix = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)
    Iy = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=3)
    if show_steps:
        plt.subplot(1, 2, 1)
        plt.imshow(Ix, cmap='gray')
        plt.title(name + ' Ix')
        plt.axis('off')
        plt.subplot(1, 2, 2)
        plt.imshow(Iy, cmap='gray')
        plt.title(name + ' Iy')
        plt.axis('off')
        plt.show()
        
    # compute Ix^2, Iy^2, Ix*Iy
    Ix2 = Ix * Ix
    Iy2 = Iy * Iy
    Ixy = Ix * Iy
    
    # W eighted average of Ix^2, Iy^2, Ix*Iy with Gaussian filter
    kernel_size = 3
    sigma = 1
    Sx2 = cv2.GaussianBlur(Ix2, (kernel_size, kernel_size), sigma)
    Sy2 = cv2.GaussianBlur(Iy2, (kernel_size, kernel_size), sigma)
    Sxy = cv2.GaussianBlur(Ixy, (kernel_size, kernel_size), sigma)
    if show_steps:
        plt.subplot(1, 3, 1)
        plt.imshow(Sx2, cmap='gray')
        plt.title(name + ' Sx2')
        plt.axis('off')
        plt.subplot(1, 3, 2)
        plt.imshow(Sy2, cmap='gray')
        plt.title(name + ' Sy2')
        plt.axis('off')
        plt.subplot(1, 3, 3)
        plt.imshow(Sxy, cmap='gray')
        plt.title(name + ' Sxy')
        plt.axis('off')
        plt.show()
    
    # Compute the Harris corner response function R
    detM = (Sx2 * Sy2) - (Sxy ** 2)
    traceM = Sx2 + Sy2
    R = detM - k * (traceM ** 2)
    if show_steps:
        plt.imshow(R, cmap='gray')
        plt.title(name + ' Response function R')
        plt.axis('off')
        plt.show()
    
    # Thresholding
    R_max = np.max(R)
    corners = np.zeros_like(gray, dtype=np.uint8)
    threshold_value = threshold * R_max
    corners[R > threshold_value] = 255
    
    # non-maximum suppression
    dst = cv2.dilate(R, None)
    corners = np.zeros_like(gray, dtype=np.uint8)
    corners[(R == dst) & (R > threshold_value)] = 255
    
    # mark the corners in the original image
    gray_rgb = cv2.cvtColor(gray, cv2.COLOR_GRAY2RGB)
    corner_points = np.argwhere(corners)
    for point in corner_points:
        y, x = point
        cv2.circle(gray_rgb, (x, y), radius=3, color=(255, 0, 0), thickness=-1) 

    if show:
        plt.imshow(gray_rgb)
        plt.title(name + ' red corner points')
        plt.axis('off')
        plt.show()
    
    return gray_rgb

corner_dark_1_ = harris_corner_detection(image_path_dark, "Dark-1", k=0.04, threshold=0.01, show_steps=True)
corner_bright_1_ = harris_corner_detection(image_path_bright, "Bright-1", k=0.04, threshold=0.01, show_steps=True)

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

[Fill in your discussion to Step 3 here]

We can observe that the detected corners in dark image result is much less than those in bright image. Besides, the corners in bright image distributes more uniformly than which in dark image. The corners of the dark image mainly concetrate on bright area, and lose the dark area's corners. However, the dark area takes most part of the picture.

Possible reason for the result of dark image is that, the gradient of the dark area are all small and hard to distinguish from this point because of low contrast. The low light also means points are easy to be surpressed by the algorithm.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 4 IS HERE                           ##
#############################################################################
import cv2
import numpy as np
import matplotlib.pyplot as plt

def histogram_equalization(image_path, output_path):
    """
    Histogram equalization for the input image.

    arguments:
    - image_path: Path to the input image file.
    - output_path: Path to save the output image file.

    return:
    - equalized_image: Image after histogram equalization.
    """
import cv2
import numpy as np

def histogram_equalization(image_path, output_path):
    """
    Histogram equalization for the input image.

    arguments:
    - image_path: Path to the input image file.
    - output_path: Path to save the output image file.

    return:
    - equalized_image: Image after histogram equalization.
    """
    image = cv2.imread(image_path)
    if image is None:
        raise FileNotFoundError("image file not found.")

    # interface for different image types
    # gray
    if len(image.shape) == 2 or image.shape[2] == 1:
        equalized_image = equalize_gray_image(image)
    # color
    elif len(image.shape) == 3 and image.shape[2] == 3:
        equalized_image = equalize_color_image(image)
    else:
        raise ValueError("unsupported image type, only support gray and color image.")
    cv2.imwrite(output_path, equalized_image)
    return equalized_image

def equalize_gray_image(gray_image):
    """
    Perform histogram equalization on a grayscale image.

    Parameters:
    - gray_image: Input grayscale image (numpy array).

    Returns:
    - equalized_gray: Grayscale image after histogram equalization.
    """
    hist, bins = np.histogram(gray_image.flatten(), 256, [0,256])

    cdf = hist.cumsum()
    cdf_normalized = cdf * (255 / cdf[-1])  # normalize
    # map the pixel value to the equalized value 
    equalized_gray = np.interp(gray_image.flatten(), bins[:-1], cdf_normalized)
    equalized_gray = equalized_gray.reshape(gray_image.shape).astype('uint8')

    return equalized_gray

def equalize_color_image(color_image):
    """
    Perform histogram equalization on a color image.

    Parameters:
    - color_image: Input color image (numpy array).
    
    Returns:
    - equalized_color: Color image after histogram equalization.
    """
    y_cr_cb = cv2.cvtColor(color_image, cv2.COLOR_BGR2YCrCb)
    y_channel, cr_channel, cb_channel = cv2.split(y_cr_cb)
    # he on brightness channel
    equalized_y = equalize_gray_image(y_channel)
    # inconvert back to BGR
    equalized_y_cr_cb = cv2.merge((equalized_y, cr_channel, cb_channel))
    equalized_color = cv2.cvtColor(equalized_y_cr_cb, cv2.COLOR_YCrCb2BGR)

    return equalized_color

if __name__ == '__main__':
    dark_he_path = 'dark-he.jpg' 
    bright_he_path = 'bright-he.jpg'

    dark_equalized_img = histogram_equalization(image_path_dark, dark_he_path)
    bright_equalized_img = histogram_equalization(image_path_bright, bright_he_path)

    # dark-he show
    original_image = cv2.imread(image_path_dark)
    plt.figure(figsize=(12, 6))
    plt.subplot(1, 2, 1)
    plt.imshow(cv2.cvtColor(original_image, cv2.COLOR_BGR2RGB))
    plt.title('Dark-1')
    plt.axis('off')

    plt.subplot(1, 2, 2)
    plt.imshow(cv2.cvtColor(dark_equalized_img, cv2.COLOR_BGR2RGB))
    plt.title('Dark-he')
    plt.axis('off')
    plt.show()
    
    # bright-he show
    original_image = cv2.imread(image_path_bright)
    plt.figure(figsize=(12, 6))
    plt.subplot(1, 2, 1)
    plt.imshow(cv2.cvtColor(original_image, cv2.COLOR_BGR2RGB))
    plt.title('Bright-1')
    plt.axis('off')

    plt.subplot(1, 2, 2)
    plt.imshow(cv2.cvtColor(bright_equalized_img, cv2.COLOR_BGR2RGB))
    plt.title('Bright-he')
    plt.axis('off')
    plt.show()

[Fill in your discussion to Step 4 here]

Histogram equalization significantly improves image contrast, making details more visible in the previously dark "Dark-1" image. By redistributing pixel intensities, the HE-enhanced image has a broader range of values, revealing textures and objects that were hidden in shadows. This improvement can make objects like plants and structures clearer.

For Canny edge detection, HE will lead to sharper, more distinct edges as the contrast between regions is enhanced. More edges will likely be detected, especially in areas that were too dark before.

Similarly, Harris corner detection will identify more corners due to the stronger intensity variations. The corners will also be more accurately localized as gradients become clearer. Overall, HE prepares the image for better feature extraction, improving both edge and corner detection outcomes.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 5 IS HERE                           ##
#############################################################################

def plot_comparison(original_image, he_image, title_original, title_he):
    plt.figure(figsize=(10, 5))

    # Original image
    plt.subplot(1, 2, 1)
    plt.imshow(cv2.cvtColor(original_image, cv2.COLOR_BGR2RGB))
    plt.title(title_original)
    plt.axis('off')

    # HE image
    plt.subplot(1, 2, 2)
    plt.imshow(cv2.cvtColor(he_image, cv2.COLOR_BGR2RGB))
    plt.title(title_he)
    plt.axis('off')

    plt.show()

corner_dark_1 = harris_corner_detection(image_path_dark, "Dark-1", k=0.04, threshold=0.01, show_steps=False, show=False)
corner_bright_1 = harris_corner_detection(image_path_bright, "Bright-1", k=0.04, threshold=0.01, show_steps=False, show=False)
edge_dark_1 = canny_edge_detection(image_path_dark, "Dark-1", low_threshold=50, high_threshold=150, show_steps=False, show=False)
edge_bright_1 = canny_edge_detection(image_path_bright, "Bright-1", low_threshold=50, high_threshold=150, show_steps=False, show=False)

edge_dark_he = canny_edge_detection(dark_he_path, "Dark-he", low_threshold=50, high_threshold=150, show_steps=False, show=False)
edge_bright_he = canny_edge_detection(bright_he_path, "Bright-he", low_threshold=50, high_threshold=150, show_steps=False, show=False)
corner_dark_he = harris_corner_detection(dark_he_path, "Dark-he", k=0.04, threshold=0.01, show_steps=False, show=False)
corner_bright_he = harris_corner_detection(bright_he_path, "Bright-he", k=0.04, threshold=0.01, show_steps=False, show=False)

# draw comparison plot
plot_comparison(edge_dark_1, edge_dark_he, 'Dark-1 Edge Detection', 'Dark-he Edge Detection')
plot_comparison(edge_bright_1, edge_bright_he, 'Bright-1 Edge Detection', 'Bright-he Edge Detection')
plot_comparison(corner_dark_1, corner_dark_he, 'Dark-1 Corner Detection', 'Dark-he Corner Detection')
plot_comparison(corner_bright_1, corner_bright_he, 'Bright-1 Corner Detection', 'Bright-he Corner Detection')

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

[Fill in your discussion to Step 5 here]

Comparison with Hypothesis:

1. Canny Edge Detection:

   Dark-1 vs. Dark-he: The initial hypothesis was that after applying HE to Dark-1 (creating Dark-he), Canny edge detection would detect more edges and those edges would appear sharper. The results confirm this, as Dark-he Edge Detection displays more visible edges compared to Dark-1 Edge Detection, which is almost completely black due to the lack of contrast. HE has clearly enhanced the contrast, making edges around objects like leaves and building structures more detectable.

2. Harris Corner Detection:

   Dark-1 vs. Dark-he: The hypothesis was that Harris corner detection would identify more corners in Dark-he, with better localization, due to the improved intensity variations after HE. The results align with this, as Dark-he Corner Detection shows a dense distribution of corners (marked in blue), while Dark-1 Corner Detection has almost little corner detection due to the lack of visible intensity differences. HE has enabled the Harris algorithm to detect more corners by enhancing the image's gradient details.

3. Rationalization of Hypothesis:

   The results support the hypothesis well, as both Canny and Harris detection significantly improved on Dark-he after histogram equalization. The increased contrast allowed both algorithms to better capture the image's features (edges and corners). For Bright-1, the results were also improved after HE, satisfying the hypothesis.

Comparison to Bright-1:

1. Canny Edge Detection: Comparing Bright-1 Edge Detection and Bright-he Edge Detection, both images show many visible edges, but there is a slight improvement in edge sharpness and clarity in Bright-he due to the histogram equalization. However, the difference is not as dramatic as it is with the Dark-1 image because Bright-1 already had enough contrast for edge detection, so HE only provided a marginal improvement.

2. Harris Corner Detection: After applying HE, Bright-he Corner Detection detects more corners, and more detailed corner in dark area, which means he also works on bright picture and improve the feature extraction.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 6 IS HERE                           ##
#############################################################################


# from PIL import Image
import matplotlib.pyplot as plt
import cv2
import numpy as np


def linear_enhance(image, alpha=1.5, beta=0):
    """
    Enhances an image using a linear transformation: new_pixel = alpha * pixel + beta.
    
    Parameters:
    - image : numpy.ndarray
        The input image to be enhanced.
    - alpha : float, optional
        The gain factor (contrast control). A value > 1 increases contrast, while a value < 1 decreases contrast.
        Default is 1.5.
    - beta : int, optional
        The bias value (brightness control). Positive values increase brightness, while negative values decrease it.
        Default is 0.
    
    Returns:
    - enhanced_image : numpy.ndarray
        The enhanced image with the same dimensions as the input image.
    """
    # Apply linear transformation to the image: new_pixel = alpha * pixel + beta
    enhanced_image = cv2.convertScaleAbs(image, alpha=alpha, beta=beta)
    
    return enhanced_image


def sharpen_image_adjustable(image, alpha=0.2):
    """
    Sharpens an image using a convolution kernel with adjustable sharpening strength.

    Parameters:
    image : numpy.ndarray
        The input image to be sharpened.
    alpha : float, optional
        The weight applied to the sharpened image. Higher values increase sharpening (default is 1.5).

    Returns:
    sharpened_image : numpy.ndarray
        The sharpened image.
    """
    # Define the sharpening kernel
    kernel = np.array([[0, 0, -1, 0, 0],
                   [0, -1, -2, -1, 0],
                   [-1, -2, 16, -2, -1],
                   [0, -1, -2, -1, 0],
                   [0, 0, -1, 0, 0]])
    
    # Apply the kernel to the image
    sharpened_image = cv2.filter2D(image, -1, kernel)
    
    # Blend the original and sharpened images using alpha for sharpening strength
    output_image = cv2.addWeighted(image, 1 - alpha, sharpened_image, alpha, 0)
    
    return output_image


def denoise_image_nl_means(image, h=10, template_window_size=7, search_window_size=21):
    """
    Denoises an image using the Non-Local Means Denoising algorithm.
    
    Parameters:
    - image : The input image to be denoised.
    - h :  Parameter regulating filter strength. A larger h means more smoothing (default is 10).
    - template_window_size : Size of the window used to compute the weights (default is 7).
    - search_window_size : Size of the window used to search for patches similar to the central pixel (default is 21).

    Returns:
    - denoised_image : The denoised image.
    """
    denoised_image = cv2.fastNlMeansDenoising(image, h=h, 
                                              templateWindowSize=template_window_size,
                                              searchWindowSize=search_window_size)
    return denoised_image

def multi_scale_retinex(img, scales=[15, 80, 250]):
    """
    Multi-Scale Retinex (MSR) algorithm to enhance the image by reducing illumination effects.
    
    Parameters:
    - img: Input image (BGR format).
    - scales: List of scales (sigma values) for Gaussian blurs in MSR. Default is [15, 80, 250].
    
    Returns:
    - msr_img: The enhanced image after MSR.
    """
    # Convert image to float32 and add 1 to avoid log(0) happens
    img = np.float32(img) + 1.0
    msr_result = np.zeros_like(img)

    for scale in scales:
        img_blur = cv2.GaussianBlur(img, (0, 0), sigmaX=scale, sigmaY=scale)
        retinex = np.log(img) - np.log(img_blur + 1.0)
        msr_result += retinex
        
    # Average the results from all scales
    msr_result /= len(scales)
    msr_normalized = cv2.normalize(msr_result, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX)
    msr_img = np.uint8(msr_normalized)

    return msr_img

def msr_enhancement(image_path, output_path, scales=[15, 80, 250]):
    """
    Apply Multi-Scale Retinex (MSR) enhancement to the input image and save the result.
    
    Parameters:
    - image_path: File path of the input image.
    - output_path: File path to save the enhanced image.
    - scales: List of scales for Gaussian blurs in MSR. Default is [15, 80, 250].
    
    Returns:
    - enhanced_img: The enhanced image after MSR.
    """
    img = cv2.imread(image_path)
    if img is None:
        print("Could not load the image. Please check the file path.")
        return

    linear_enhanced_image = linear_enhance(img, alpha=1.8, beta=20)
    # Apply MSR to each channel (R, G, B)
    enhanced_channels = []
    for i in range(3):
        enhanced_channel = multi_scale_retinex(linear_enhanced_image[:, :, i], scales=scales)
        enhanced_channels.append(enhanced_channel)
    enhanced_img = cv2.merge(enhanced_channels)
    # Sharpen the denoised image
    enhanced_img = sharpen_image_adjustable(enhanced_img, alpha=0.2)
    # Denoise the enhanced image
    enhanced_img = denoise_image_nl_means(enhanced_img, h=7)
    cv2.imwrite(output_path, enhanced_img)

    return enhanced_img


retinex_image_path = 'dark-retinex.jpg'
# Apply mix enhancement (including linear-enhancement, multi-scale-retinex, sharpening and denoising) to the dark image
enhanced_img = msr_enhancement(image_path_dark, retinex_image_path, scales=[70, 200])

if enhanced_img is not None:
    img_rgb = cv2.cvtColor(enhanced_img, cv2.COLOR_BGR2RGB)
    plt.imshow(img_rgb)
    plt.title("Dark-1: Multi-Scale Retinex Enhanced Image")
    plt.axis('off')
    plt.show()

edge_dark_retinex = canny_edge_detection(retinex_image_path, "Dark-retinex edge detection", low_threshold=50, high_threshold=150, show_steps=False, show=False)
corner_dark_retinex = harris_corner_detection(retinex_image_path, "Dark-retinex corner detection", k=0.04, threshold=0.01, show_steps=False, show=False)

# draw comparison plot
plot_comparison(edge_dark_retinex, edge_dark_1, 'Dark-self Edge Detection', 'Dark-1 Edge Detection')
plot_comparison(edge_dark_retinex, edge_dark_he, 'Dark-self Edge Detection', 'Dark-he Edge Detection')
plot_comparison(edge_dark_retinex, edge_bright_1, 'Dark-self Edge Detection', 'Bright-1 Edge Detection')
plot_comparison(corner_dark_retinex, corner_dark_1, 'Dark-self Corner Detection', 'Dark-1 Corner Detection')
plot_comparison(corner_dark_retinex, corner_dark_he, 'Dark-self Corner Detection', 'Dark-he Corner Detection')
plot_comparison(corner_dark_retinex, corner_bright_he, 'Dark-self Corner Detection', 'Bright-he Corner Detection')

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

[Fill in your discussion to Step 6 here]

I applied the multi-retinex enhancement algorithm, which try to remove reflection from the picture and show the origin brightness of objects.

Comparing Dark-1, Dark-1-HE, Dark-1-Self, and Bright-1:

• Edge Detection:

  • Dark-1 shows minimal edge detection due to its low contrast. Dark-1-HE significantly improves edge detection by enhancing contrast but can lead to unnatural sharp edges, especially in previously dark areas. Dark-1-Self, using multi-retinex, offers a more balanced enhancement, revealing edges in dark regions without over-exaggerating contrasts. Bright-HE already had good edge detection due to its inherent brightness, and both Bright-HE and Dark-1-Self exhibit smoother transitions.
  • Dark-1-Self, enhanced by multi-retinex, captures a more natural lighting balance than Dark-1-HE, where histogram equalization can sometimes over-enhance darker areas, resulting in unnatural contrast. The edges detected in Dark-1-Self are smoother and more reflective of real-world lighting conditions, while Dark-1-HE provides sharper, but sometimes overly enhanced, edges. The multi-retinex method excels in preserving subtle textures while revealing hidden details in shadowed regions.

• Corner Detection:

  • Dark-1 detects very few corners because of its low contrast. After HE, Dark-1-HE improves corner detection by highlighting strong intensity changes, but it can exaggerate areas of noise. Dark-1-Self, with the multi-retinex method, provides a more natural distribution of corners, better capturing fine details in both shadowed and bright regions. It maintains a balance, avoiding excessive noise while revealing useful structural information. Bright-1 performs well in corner detection naturally due to good initial lighting, though Dark-1-Self shows comparable performance in bringing out details from the originally darker regions.
  • In Dark-1-Self, the corners detected are more dispersed. The multi-retinex enhancement maintains a balance between dark and bright regions, resulting in a more detailed corner map. Dark-1-HE, in contrast, tends to boost corners in areas with high-intensity changes, sometimes overemphasizing areas where noise is introduced through excessive contrast adjustments.

• Rationalization:

  1. Dark-1: Low contrast results in minimal edge and corner detection, with most features hidden due to darkness.

  2. Dark-1-HE: Histogram equalization boosts contrast significantly, revealing more edges and corners, but often leads to overly sharp, exaggerated features and noise.

  3. Dark-1-Self (Multi-Retinex): Provides a more balanced enhancement compared to HE, revealing edges and corners naturally without over-exaggeration. It provides smooth transitions, avoiding noise while still enhancing darker regions.

  4. Bright-1: Naturally performs well in both edge and corner detection due to good contrast and lighting.

Before Task 2

• For the following task, you are to use GT1 (denote as Bright-1) and GT2 (denote as Bright-2) as the stereo image pair.
• All raw images are given in RGB format (colored images), but you may also convert to grayscale images for your convenience. Show the step and result of grayscale conversion first if you are to convert to grayscale for your tasks.
• You may opt to resize the images in case you think the image is too large, the resize images must not be smaller than 800*800. Show the step and result of resize first if you opt to resize the images.
• IMPORTANT! You may use any function of OpenCV or its equivalence for basic operations (e.g., loading images, matrix computation, etc.), but strictly NOT the direct functions for each individual task/step (e.g., cv.computeCorrespondEpilines for computing and drawing the epipolar lines and cv.findFundamentalMat for computing the fundamental matrix). Using such functions would consider the code to be erroneous.

Task 2: Computing the Fundamental Matrix and Finding the Epipolar Lines for Stereo Image Pairs (40%)

In this task, you will need to compute the fundamental matrix based on the (Normalized) 8-point algorithm. You are to choose the keypoints manually (you can use tools such as Paint for Windows), and then compute the fundamental matrix. You would then find, draw out, observe and discuss about the epipolar lines for your chosen keypoints. You are to follow the following steps (requirements):

Steps/Requirements for Task 2:

1. Find and highlight the keypoints in the image manually. For the 8-point algorithm, you should select at least 8 non-colinear points. State the individual positions in homogeneous coordinates, and highlight them in the corresponding image. If you choose to use the normalized 8-pt algorithm, you should state the coordinates of both the original keypoints and the normalized keypoints and showcase your code for normalization, but you only need to highlight the original keypoints. If you find highlighting in the pixel level is difficult, you may use a small circle centered at the selected keypoint for highlighting. (10%)
2. With the selected keypoints, compute the fundamental matrix F in detail. State the obtained fundamental matrix. (10%)
3. With the fundamental matrix, draw all the epipolar lines corresponding to the selected keypoints on both images. Observe and discuss how these lines indicate the positioning of the cameras. Also discuss if the computed lines match your intuitive idea of how the lines should be formed. If yes, why? If no, why not? (10%)
4. Lastly, with the computed fundamental matrix, we would be also able to find the epipolar line and subsequently the corresponding keypoint given a new keypoint. Select a new keypoint that does not lie on any of the drawn epipolar lines on Bright-1, then draw the epipolar line corresponding to this new keypoint on Bright-2. Observe and discuss if a possible keypoint can be obtained by searching across the drawn epipolar line. (10%)

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 1 IS HERE                           ##
#############################################################################

import cv2
import json
from functools import partial
import matplotlib.pyplot as plt

# 回调函数，用于处理鼠标点击事件
def click_event(event, x, y, flags, param):
    img, container, img_id, click_count = param
    if event == cv2.EVENT_LBUTTONDOWN:
        click_count[0] += 1
        print(f"Clicked at: ({x}, {y})")
        container.append((x, y))
        font = cv2.FONT_HERSHEY_SIMPLEX
        cv2.putText(img, f"({x}, {y})  #{click_count[0]}", (x, y), font, 0.5, (255, 0, 0), 1, cv2.LINE_AA)
        cv2.imshow('GT1 Image', img)
        if len(container) >= 8:
            with open(f'GT{img_id}_coordinates.json', 'w') as f:
                json.dump(container, f)
            cv2.putText(img, "dump json", (400, 400), font, 0.5, (255, 0, 0), 1, cv2.LINE_AA)

# 读取图像
GT1_img = cv2.imread('GT1.jpg')
GT1_img = cv2.resize(GT1_img, (800, 800))
click_count = [0]
GT1_coordinates = []

# 显示图像
plt.imshow(cv2.cvtColor(GT1_img, cv2.COLOR_BGR2RGB))
plt.axis('off')
plt.title("GT1 Image")


# Display the images using OpenCV
cv2.imshow('GT1 Image', GT1_img)
cv2.setMouseCallback('GT1 Image', click_event, (GT1_img, GT1_coordinates, 1, click_count))
cv2.waitKey(0)
cv2.destroyAllWindows()

#############################################################################
##                                                                         ##
##                    TODO: CONTINUATION FOR STEP 1                        ##
#############################################################################

import cv2
import json
from functools import partial
import matplotlib.pyplot as plt

def click_event(event, x, y, flags, param):
    img, container, img_id, click_count = param
    if event == cv2.EVENT_LBUTTONDOWN:
        click_count[0] += 1
        print(f"Clicked at: ({x}, {y})")
        container.append((x, y))
        font = cv2.FONT_HERSHEY_SIMPLEX
        cv2.putText(img, f"({x}, {y})  #{click_count[0]}", (x, y), font, 0.5, (255, 0, 0), 1, cv2.LINE_AA)
        cv2.imshow('GT2 Image', img)
        if len(container) >= 8:
            with open(f'GT{img_id}_coordinates.json', 'w') as f:
                json.dump(container, f)
            cv2.putText(img, "dump json", (400, 400), font, 0.5, (255, 0, 0), 1, cv2.LINE_AA)

GT2_img = cv2.imread('GT2.jpg')
GT2_img = cv2.resize(GT2_img, (800, 800))
click_count = [0]
GT2_coordinates = []

plt.imshow(cv2.cvtColor(GT2_img, cv2.COLOR_BGR2RGB))
plt.axis('off')
plt.title("GT2 Image")


# Display the images using OpenCV
cv2.imshow('GT2 Image', GT2_img)
cv2.setMouseCallback('GT2 Image', click_event, (GT2_img, GT2_coordinates, 2, click_count))
cv2.waitKey(0)
cv2.destroyAllWindows()

#############################################################################
##                                                                         ##
##                    TODO: CONTINUATION FOR STEP 1                        ##
#############################################################################

import cv2
import matplotlib.pyplot as plt
import json
import numpy as np

def to_homogeneous(coords):
    """Convert coordinates to homogeneous coordinates."""
    return np.hstack([coords, np.ones((coords.shape[0], 1))])

def normalize_points(points):
    """Normalize the points."""
    mean = np.mean(points, axis=0)
    std = np.std(points)
    T = np.array([[1/std, 0, -mean[0]/std],
                  [0, 1/std, -mean[1]/std],
                  [0, 0, 1]])
    points_normalized = np.dot(T, points.T).T
    return points_normalized, T

# 读取图像
image_1 = cv2.imread('image_select_gt1.png')
image_1_rgb = cv2.cvtColor(image_1, cv2.COLOR_BGR2RGB)

image_2 = cv2.imread('image_select_gt2.png')
image_2_rgb = cv2.cvtColor(image_2, cv2.COLOR_BGR2RGB)

# 读取 JSON 文件中的坐标
with open('GT1_coordinates.json', 'r') as f:
    gt1_coordinates = np.array(json.load(f))

with open('GT2_coordinates.json', 'r') as f:
    gt2_coordinates = np.array(json.load(f))

# 将坐标转换为齐次坐标
gt1_homogeneous = to_homogeneous(gt1_coordinates)
gt2_homogeneous = to_homogeneous(gt2_coordinates)

# 归一化坐标
gt1_normalized, T1 = normalize_points(gt1_homogeneous)
gt2_normalized, T2 = normalize_points(gt2_homogeneous)

# 在图像上绘制小红圈
for (x, y) in gt1_coordinates:
    cv2.circle(image_1_rgb, (int(x), int(y)), 5, (255, 0, 0), -1)

for (x, y) in gt2_coordinates:
    cv2.circle(image_2_rgb, (int(x), int(y)), 5, (255, 0, 0), -1)

# 显示图像
plt.subplot(1, 2, 1)
plt.imshow(image_1_rgb)
plt.axis('off')
plt.title('Image GT1 with Keypoints')

plt.subplot(1, 2, 2)
plt.imshow(image_2_rgb)
plt.axis('off')
plt.title('Image GT2 with Keypoints')
plt.show()

# 打印原始和归一化后的坐标
print("Original GT1 Coordinates (Homogeneous):")
print(gt1_homogeneous)
print("Normalized GT1 Coordinates:")
print(gt1_normalized)

print("Original GT2 Coordinates (Homogeneous):")
print(gt2_homogeneous)
print("Normalized GT2 Coordinates:")
print(gt2_normalized)

Original GT1 Coordinates (Homogeneous):
[[265.  40.   1.]
 [115. 126.   1.]
 [260. 169.   1.]
 [584. 303.   1.]
 [438. 341.   1.]
 [647. 570.   1.]
 [232. 634.   1.]
 [598. 640.   1.]
 [464. 385.   1.]]
Normalized GT1 Coordinates:
[[-0.5656706  -1.32268461  1.        ]
 [-1.1926454  -0.96321906  1.        ]
 [-0.58656976 -0.78348629  1.        ]
 [ 0.76769581 -0.2233888   1.        ]
 [ 0.15744034 -0.06455518  1.        ]
 [ 1.03102523  0.89262634  1.        ]
 [-0.70360505  1.16013559  1.        ]
 [ 0.82621346  1.18521458  1.        ]
 [ 0.26611597  0.11935742  1.        ]]
Original GT2 Coordinates (Homogeneous):
[[364.  45.   1.]
 [199. 136.   1.]
 [356. 170.   1.]
 [553. 307.   1.]
 [390. 343.   1.]
 [390. 586.   1.]
 [147. 610.   1.]
 [414. 658.   1.]
 [410. 389.   1.]]
Normalized GT2 Coordinates:
[[ 0.02706998 -1.45003172  1.        ]
 [-0.73140029 -1.03172387  1.        ]
 [-0.00970433 -0.87543303  1.        ]
 [ 0.8958632  -0.24567286  1.        ]
 [ 0.14658651 -0.08018844  1.        ]
 [ 0.14658651  1.03683142  1.        ]
 [-0.97043335  1.14715436  1.        ]
 [ 0.25690946  1.36780026  1.        ]
 [ 0.2385223   0.13126388  1.        ]]

[Fill in your discussion to Step 1 here]

The corresponding points and coordinates are shown on picture and printed.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 2 IS HERE                           ##
#############################################################################

import cv2
import numpy as np
import json

def to_homogeneous(coords):
    """Convert coordinates to homogeneous coordinates."""
    return np.hstack([coords, np.ones((coords.shape[0], 1))])

def normalize_points(points):
    """Normalize the points."""
    mean = np.mean(points, axis=0)
    std = np.std(points)
    T = np.array([[1/std, 0, -mean[0]/std],
                  [0, 1/std, -mean[1]/std],
                  [0, 0, 1]])
    points_normalized = np.dot(T, points.T).T
    return points_normalized, T

def construct_matrix_A(points1, points2):
    """Construct matrix A for the 8-point algorithm."""
    A = []
    for i in range(points1.shape[0]):
        x1, y1 = points1[i, 0], points1[i, 1]
        x2, y2 = points2[i, 0], points2[i, 1]
        A.append([x2*x1, x2*y1, x2, y2*x1, y2*y1, y2, x1, y1, 1])
    return np.array(A)

def compute_fundamental_matrix(points1, points2):
    """Compute the fundamental matrix using the normalized 8-point algorithm."""
    points1_homogeneous = to_homogeneous(points1)
    points2_homogeneous = to_homogeneous(points2)
    
    points1_normalized, T1 = normalize_points(points1_homogeneous)
    points2_normalized, T2 = normalize_points(points2_homogeneous)
    
    A = construct_matrix_A(points1_normalized, points2_normalized)
    
    U, S, Vt = np.linalg.svd(A)
    F_normalized = Vt[-1].reshape(3, 3)
    
    U, S, Vt = np.linalg.svd(F_normalized)
    S[2] = 0
    F_normalized = np.dot(U, np.dot(np.diag(S), Vt))
    
    F = np.dot(T2.T, np.dot(F_normalized, T1))
    
    return F / F[2, 2]

with open('GT1_coordinates.json', 'r') as f:
    gt1_coordinates = np.array(json.load(f))
gt1_coordinates = gt1_coordinates[:7]+gt1_coordinates[8:]

with open('GT2_coordinates.json', 'r') as f:
    gt2_coordinates = np.array(json.load(f))
gt2_coordinates = gt2_coordinates[:7]+gt2_coordinates[8:]


# compute the fundamental matrix
F = compute_fundamental_matrix(gt1_coordinates, gt2_coordinates)

# print the fundamental matrix
print("Fundamental Matrix F:")
print(F)

Fundamental Matrix F:
[[-5.48309060e-07  3.77767136e-06 -1.41928292e-03]
 [-1.66349542e-06  1.54320564e-06  3.82817088e-04]
 [ 1.05008292e-03 -3.60053204e-03  1.00000000e+00]]

[Fill in your discussion to Step 2 here]

The corrsponding matrix is printed above.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 3 IS HERE                           ##
#############################################################################

import cv2
import numpy as np
import matplotlib.pyplot as plt
import json

def compute_correspond_epilines(points, image_id, F):
    """
    Compute the corresponding epilines for the given points.
    
    Args:
    - points (np.ndarray): Array of points in the form (N, 2).
    - which_image (int): 1 if points are from the first image, 2 if points are from the second image.
    - F (np.ndarray): Fundamental matrix.
    
    Returns:
    - epilines (np.ndarray): Array of epilines in the form (N, 3).
    """
    # Convert points to homogeneous coordinates
    points_homogeneous = np.hstack([points, np.ones((points.shape[0], 1))])
    
    if image_id == 1:
        # Compute epilines in the second image
        epilines = np.dot(F, points_homogeneous.T).T
    elif image_id == 2:
        # Compute epilines in the first image
        epilines = np.dot(F.T, points_homogeneous.T).T
    else:
        raise ValueError("which_image should be 1 or 2.")
    
    # Normalize the epilines
    epilines /= np.sqrt(epilines[:, 0]**2 + epilines[:, 1]**2).reshape(-1, 1)
    
    return epilines

def draw_lines(img, lines, points):
    r, c, _ = img.shape
    for r, pt in zip(lines, points):
        color = tuple(np.random.randint(0, 255, 3).tolist())
        x0, y0 = map(int, [0, -r[2] / r[1]])
        x1, y1 = map(int, [c, -(r[2] + r[0] * c) / r[1]])
        img = cv2.line(img, (x0, y0), (x1, y1), color, 1)
        img = cv2.circle(img, tuple(pt), 5, color, -1)
    return img

def draw_epipolar_lines(img1, img2, points1, points2, F):
    """Draw epipolar lines on the images."""
    
    # Find epilines corresponding to points in the second image (lines in the first image)
    lines1 = compute_correspond_epilines(points2, 1, F)
    img1_with_lines = draw_lines(img1, lines1, points1)

    # Find epilines corresponding to points in the first image (lines in the second image)
    lines2 = compute_correspond_epilines(points1, 2, F)
    img2_with_lines = draw_lines(img2, lines2, points2)

    return img1_with_lines, img2_with_lines

# Read the images
image_1 = cv2.imread('image_select_gt1.png')
image_1_rgb = cv2.cvtColor(image_1, cv2.COLOR_BGR2RGB)

image_2 = cv2.imread('image_select_gt2.png')
image_2_rgb = cv2.cvtColor(image_2, cv2.COLOR_BGR2RGB)

# Read the coordinates
with open('GT1_coordinates.json', 'r') as f:
    gt1_coordinates = np.array(json.load(f))

with open('GT2_coordinates.json', 'r') as f:
    gt2_coordinates = np.array(json.load(f))

img1_with_lines, img2_with_lines = draw_epipolar_lines(image_1_rgb, image_2_rgb, gt1_coordinates, gt2_coordinates, F)
print(f"Fundamental Matrix F: \n{F}")

# Display the images with epipolar lines
plt.figure(figsize=(15, 10))

plt.subplot(1, 2, 1)
plt.imshow(img1_with_lines)
plt.axis('off')
plt.title('Image GT1 with Epipolar Lines')

plt.subplot(1, 2, 2)
plt.imshow(img2_with_lines)
plt.axis('off')
plt.title('Image GT2 with Epipolar Lines')
plt.show()

Fundamental Matrix F: 
[[-5.48309060e-07  3.77767136e-06 -1.41928292e-03]
 [-1.66349542e-06  1.54320564e-06  3.82817088e-04]
 [ 1.05008292e-03 -3.60053204e-03  1.00000000e+00]]

[Fill in your discussion to Step 3 here]

The epipolar lines should ideally radiate from the epipole and pass through the corresponding points in both images. However, The provided image shows epipolar lines that do not perfectly pass through the epipoles. This mismatch between expectation and result occurs due to several possible reasons:

1. Inaccuracy in Point Matching:

   • The accuracy of the point correspondences between the two images plays a crucial role in determining the correctness of the epipolar geometry. If the matching points are noisy or inaccurately placed, the computed fundamental matrix (or essential matrix) might not reflect the true epipolar geometry, causing the epipolar lines to deviate from the epipoles.

2. Imperfect Estimation of the Fundamental Matrix:

   • The fundamental matrix is estimated based on the corresponding points between the two images. If these points are not optimally chosen or if there are outliers, the resulting matrix may not be accurate. This would lead to incorrect epipolar lines that fail to converge at the correct epipole.

3. Numerical Instability and Precision:

   • The algorithm used (such as the eight-point algorithm) may suffer from numerical instability, especially if the point coordinates are not well-matched between the pair of stero image or properly alligned. This can lead to a suboptimal estimation of the fundamental matrix and cause the epipolar lines to deviate from where they should ideally be.

#############################################################################
##                                                                         ##
##           TODO: CODE BLOCK FOR STEP 4 IS HERE                           ##
#############################################################################

import cv2
import numpy as np
import matplotlib.pyplot as plt
import json

def draw_lines_pure(img, lines):
    r, c, _ = img.shape
    for r in lines:
        color = tuple(np.random.randint(0, 255, 3).tolist())
        x0, y0 = map(int, [0, -r[2] / r[1]])
        x1, y1 = map(int, [c, -(r[2] + r[0] * c) / r[1]])
        img = cv2.line(img, (x0, y0), (x1, y1), color, 1)
    return img

def draw_circle(img, points):
    for point in points:
        img = cv2.circle(img, tuple(point), 5, (255, 0, 0), -1)
    return img

image_2 = cv2.imread('GT2.jpg')
image_2 = cv2.resize(image_2, (800, 800))
image_2_rgb = cv2.cvtColor(image_2, cv2.COLOR_BGR2RGB)

new_keypoint = np.array([[200, 300]])
img1_with_new_point = draw_circle(img1_with_lines.copy(), new_keypoint)
new_epiline = compute_correspond_epilines(new_keypoint, 1, F)
img2_with_new_epiline = draw_lines_pure(image_2_rgb.copy(), new_epiline)

plt.figure(figsize=(15, 10))

plt.subplot(1, 2, 1)
plt.imshow(img1_with_new_point)
plt.axis('on')
plt.title('Image GT1 with Epipolar Lines')

plt.subplot(1, 2, 2)
plt.imshow(img2_with_new_epiline)
plt.axis('on')
plt.title('Image GT2 with New Epipolar Line')
plt.show()

[Fill in your discussion to Step 4 here]

It seems that due to the inaccuracy in previous calculated suboptimal fundemental matrix, the drawn corresponding epipolar is not as expected accurate, make it hard to find solution on the corresponding epipolar line by using epipolar constraint and search on the epipolar line.